# Libraries
library(tidyverse)
library(lubridate)
library(corrr)
library(corrplot)
library(ggpubr)
library(modelr)
library(effectsize)
library(glmpath)
library(glmnet)
library(MASS)
library(ggcorrplot)
library(broom)
library(ggfortify)
library(relaimpo)
library(olsrr)
library(parameters)
library(BayesFactor)
# for meff function
source('https://osf.io/ws6xc/download')
# Files
home <- "~/Box/Mooddata_Coordinating/BABIES/Data/final_scored_data/"
eligible_id_file <- "~/Box/lucy_king_files/BABIES/lena_symptoms/inclusion_tracker_sub2.xlsx"
lena_edited_file <- "~/Box/lucy_king_files/BABIES/lena_symptoms/included_recording_dates.csv"
fp_intervention_id <- paste0(home, "lab_caregiving_behavior/free_play_intervention_assignment.csv")
itsea_file <- paste0(home, "survey_18month/itsea_scored_20201001.csv")
demographics_file <- paste0(home, "demographics/demo_6mo_cleaned_final.csv")
lena_file <- paste0(home, "LENA/lena_final_wf_day1_20200208_upto9months.csv")
ibq_file <- paste0(home, "IBQ/ibq_scored_final.csv")
sfp_care_file <- paste0(home, "lab_caregiving_behavior/PCIRS_sfp_wf_complete.csv")
fp_care_file <- paste0(home, "lab_caregiving_behavior/free_play_wf_8min_20200507.csv")
cesd_file <- paste0(home, "CESD/cesd_wf_20201001.csv")
# Functions
source("identify_outliers_histogram.R")
theme_lena <-
theme_minimal() +
theme(
panel.grid = element_blank(),
plot.title = element_text(size = 18, hjust = .5),
axis.title = element_text(size = 18),
axis.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.text = element_text(size = 16),
legend.position = "bottom"
)
# Display options
options(scipen = 999)
d0 <-
readxl::read_xlsx(eligible_id_file, sheet = "data_filtered_T3") %>%
dplyr::select(
ID,
eligible,
included
) %>%
filter(included == 1) %>%
left_join(
read_csv(lena_edited_file) %>%
dplyr::select(
ID, LENA_edited
),
by = "ID"
) %>%
left_join(
read_csv(demographics_file) %>%
dplyr::select(
ID,
baby_dob,
education,
education_txt,
annual_income_txt,
income_needs,
male,
mom_age,
mom_latinx,
mom_race,
age_behav,
secondlang,
secondlang_type
),
by = "ID"
) %>%
left_join(
read_csv(itsea_file) %>%
dplyr::select(
ID,
survey_date,
itsea_symptoms,
itsea_activity,
itsea_agress,
itsea_gad,
itsea_negemo,
itsea_depp,
itsea_intern,
itsea_extern,
# itsea symptoms items
itsea_a2,
itsea_a4,
itsea_b8,
itsea_b10,
itsea_b27,
itsea_b32,
itsea_a28,
itsea_a30,
itsea_a33,
itsea_b3,
itsea_b9,
itsea_b16,
itsea_b30,
itsea_b34,
itsea_b44,
itsea_b88,
itsea_c2,
itsea_e4,
itsea_a9,
itsea_a39,
itsea_a40,
itsea_b43,
itsea_b76,
itsea_b81,
itsea_b84,
itsea_b91,
itsea_b92,
itsea_a3,
itsea_a12,
itsea_a35,
itsea_b37,
itsea_b70,
itsea_b86,
itsea_e5,
itsea_e6,
itsea_e11,
itsea_a7,
itsea_a21,
itsea_a23,
itsea_b31,
itsea_b45,
itsea_b50,
itsea_b53,
itsea_b59,
itsea_b65,
itsea_b66,
itsea_b74,
itsea_b80,
itsea_b85,
itsea_concern_ied
),
by = "ID"
) %>%
mutate(
survey_date = parse_date_time(survey_date, orders = c("mdy HM")),
age_18mo = (baby_dob %--% survey_date) / months(1)
) %>%
left_join(
read_csv(fp_care_file) %>%
dplyr::select(ID, negmood_FP, sens_FP),
by = "ID"
) %>%
left_join(
read_csv(sfp_care_file) %>%
dplyr::select(ID, sens_R_M, sens_M),
by = "ID"
) %>%
left_join(
read_csv(ibq_file) %>%
dplyr::select(
ID,
#IBQ NEG items
ibq_64, ibq_74, ibq_75, ibq_32, ibq_79, ibq_80,
ibq_2, ibq_3_r, ibq_4, ibq_21, ibq_52, ibq_53, ibq_62,
ibq_22, ibq_76, ibq_77, ibq_78, ibq_87, ibq_89,
ibq_36, ibq_37, ibq_38, ibq_63, ibq_71, ibq_72,
NEG
),
by = "ID"
) %>%
left_join(read_csv(lena_file), by = "ID") %>%
left_join(read_csv(cesd_file), by = "ID") %>%
left_join(read_csv(fp_intervention_id), by = "ID") %>%
rename(
age_lena_d1 = age,
sens_r_sfp = sens_R_M,
sens_sfp = sens_M,
sens_fp = sens_FP,
negmood_fp = negmood_FP
)
Parsed with column specification:
cols(
ID = col_double(),
date_record = col_character(),
LENA_edited = col_double()
)
Parsed with column specification:
cols(
.default = col_double(),
due_date = col_datetime(format = ""),
baby_dob = col_datetime(format = ""),
mom_dob = col_datetime(format = ""),
baby_race = col_character(),
baby_race_describe = col_character(),
mom_race = col_character(),
momrace_describe = col_character(),
annual_income_txt = col_character(),
education_txt = col_character(),
employment_status_txt = col_character(),
marital_status_txt = col_character(),
partner_educ_txt = col_character(),
partner_employ_txt = col_character(),
questionnaire_only_date = col_date(format = ""),
educ_describe = col_logical(),
employment_explain = col_character(),
mom_pob = col_character(),
mom_native_lang = col_character(),
primarylang = col_character(),
secondlang_type = col_character()
# ... with 24 more columns
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
.default = col_double(),
survey_date = col_character(),
baby_dob = col_datetime(format = ""),
activity_concern = col_character(),
agress_concern = col_character(),
depp_concern = col_character(),
gad_concern = col_character(),
negemo_concern = col_character(),
play_concern = col_character(),
empathy_concern = col_character(),
social_concern = col_character(),
itsea_concern = col_character(),
itsea_concern_ied = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
ID = col_double(),
sens_FP = col_double(),
intrus_FP = col_double(),
posreg_FP = col_double(),
stim_FP = col_double(),
negreg_FP = col_double(),
negmood_FP = col_double(),
posmood_FP = col_double(),
detach_FP = col_double()
)
Parsed with column specification:
cols(
.default = col_double()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
.default = col_double(),
ibq_timestamp = col_datetime(format = "")
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
.default = col_double(),
date_record = col_date(format = ""),
day_type = col_character(),
lena_recordno = col_character(),
lena_probyes = col_character(),
other_caregiver = col_character(),
lena_notes = col_character()
)
See spec(...) for full column specifications.
Parsed with column specification:
cols(
ID = col_double(),
cesd_t1 = col_double(),
cesd_t2 = col_double(),
cesd_t3 = col_double(),
cesd_t4 = col_double(),
cesd_t1_concern = col_character(),
cesd_t2_concern = col_character(),
cesd_t3_concern = col_character(),
cesd_t4_concern = col_character()
)
Parsed with column specification:
cols(
ID = col_double(),
group = col_character()
)
Dyads who completed the 6-month assessment (either lab or home observation), were eligible for the 18-month assessment (i.e., were age 18 months by the time of analysis), and completed the 18-month assessment are included in the current analyses.
#missing SFP
d0 %>%
count(is.na(sens_sfp))
#missing FP
d0 %>%
count(is.na(sens_fp))
#missing IBQ
d0 %>%
count(is.na(NEG))
d0 %>%
dplyr::select(
awc_prop,
ctc_prop,
awc_hour_max,
ctc_hour_max,
sens_r_sfp,
sens_fp,
NEG,
) %>%
pair_n()
awc_prop ctc_prop awc_hour_max ctc_hour_max sens_r_sfp sens_fp NEG
awc_prop 100 100 100 100 94 96 96
ctc_prop 100 100 100 100 94 96 96
awc_hour_max 100 100 100 100 94 96 96
ctc_hour_max 100 100 100 100 94 96 96
sens_r_sfp 94 94 94 94 94 93 90
sens_fp 96 96 96 96 93 96 92
NEG 96 96 96 96 90 92 96
attr(,"class")
[1] "n_mat" "matrix"
d0 %>%
dplyr::select(ibq_64:ibq_72) %>%
psych::alpha(check.keys = TRUE)
Some items were negatively correlated with total scale and were automatically reversed.
This is indicated by a negative sign for the variable name.
Reliability analysis
Call: psych::alpha(x = ., check.keys = TRUE)
lower alpha upper 95% confidence boundaries
0.84 0.87 0.91
Reliability if an item is dropped:
Item statistics
Non missing response frequency for each item
1 2 3 4 5 6 7 miss
ibq_64 0.07 0.20 0.16 0.12 0.12 0.20 0.13 0.10
ibq_74 0.04 0.12 0.13 0.16 0.23 0.24 0.07 0.06
ibq_75 0.29 0.35 0.11 0.06 0.10 0.08 0.01 0.11
ibq_32 0.44 0.45 0.06 0.03 0.00 0.01 0.00 0.07
ibq_79 0.08 0.21 0.13 0.28 0.09 0.17 0.04 0.08
ibq_80 0.03 0.24 0.20 0.23 0.11 0.18 0.01 0.07
ibq_2 0.05 0.17 0.20 0.27 0.12 0.16 0.03 0.05
ibq_3_r 0.01 0.18 0.19 0.32 0.09 0.18 0.03 0.05
ibq_4 0.00 0.22 0.17 0.22 0.16 0.18 0.06 0.04
ibq_21 0.02 0.23 0.25 0.19 0.16 0.12 0.03 0.04
ibq_52 0.02 0.09 0.22 0.18 0.17 0.26 0.05 0.05
ibq_53 0.34 0.27 0.14 0.06 0.12 0.04 0.03 0.06
ibq_62 0.04 0.22 0.14 0.21 0.11 0.20 0.07 0.10
ibq_22 0.07 0.25 0.11 0.17 0.17 0.18 0.06 0.11
ibq_76 0.24 0.30 0.17 0.08 0.05 0.10 0.06 0.07
ibq_77 0.47 0.26 0.06 0.06 0.05 0.08 0.03 0.12
ibq_78 0.51 0.30 0.08 0.06 0.02 0.02 0.01 0.11
ibq_87 0.61 0.26 0.09 0.00 0.04 0.00 0.00 0.11
ibq_89 0.47 0.31 0.07 0.03 0.01 0.08 0.03 0.28
ibq_36 0.00 0.07 0.23 0.15 0.18 0.33 0.04 0.05
ibq_37 0.01 0.40 0.29 0.16 0.06 0.04 0.03 0.05
ibq_38 0.02 0.04 0.12 0.17 0.18 0.45 0.02 0.04
ibq_63 0.01 0.09 0.11 0.19 0.22 0.37 0.01 0.05
ibq_71 0.00 0.02 0.01 0.11 0.16 0.42 0.28 0.07
ibq_72 0.55 0.37 0.05 0.00 0.01 0.01 0.00 0.06
d0 %>%
dplyr::select(itsea_a2:itsea_b85, -itsea_c2) %>%
psych::alpha()
Item = itsea_e4 had no variance and was deleted but still is counted in the scoreSome items were negatively correlated with the total scale and probably
should be reversed.
To do this, run the function again with the 'check.keys=TRUE' option
Some items ( itsea_b84 itsea_a35 itsea_b37 itsea_e5 ) were negatively correlated with the total scale and
probably should be reversed.
To do this, run the function again with the 'check.keys=TRUE' option
Reliability analysis
Call: psych::alpha(x = .)
lower alpha upper 95% confidence boundaries
0.79 0.83 0.88
Reliability if an item is dropped:
Item statistics
Non missing response frequency for each item
0 1 2 miss
itsea_a2 0.70 0.29 0.01 0.00
itsea_a4 0.46 0.41 0.13 0.00
itsea_b8 0.31 0.51 0.18 0.00
itsea_b10 0.19 0.56 0.25 0.00
itsea_b27 0.55 0.39 0.06 0.00
itsea_b32 0.82 0.14 0.04 0.01
itsea_a28 0.80 0.19 0.01 0.00
itsea_a30 0.72 0.28 0.00 0.00
itsea_a33 0.90 0.09 0.01 0.00
itsea_b3 0.53 0.45 0.02 0.00
itsea_b9 0.64 0.32 0.04 0.00
itsea_b16 0.67 0.31 0.02 0.00
itsea_b30 0.47 0.51 0.02 0.00
itsea_b34 0.92 0.07 0.01 0.00
itsea_b44 0.19 0.63 0.18 0.00
itsea_b88 0.41 0.54 0.05 0.00
itsea_a9 0.97 0.03 0.00 0.00
itsea_a39 0.94 0.06 0.00 0.00
itsea_a40 0.98 0.02 0.00 0.00
itsea_b43 0.78 0.20 0.02 0.01
itsea_b76 0.96 0.04 0.00 0.02
itsea_b81 0.97 0.03 0.00 0.00
itsea_b84 0.98 0.01 0.01 0.00
itsea_b91 0.99 0.01 0.00 0.00
itsea_b92 0.98 0.02 0.00 0.00
itsea_a3 0.89 0.10 0.01 0.00
itsea_a12 0.83 0.15 0.02 0.01
itsea_a35 0.80 0.17 0.03 0.00
itsea_b37 0.98 0.02 0.00 0.00
itsea_b70 0.78 0.19 0.03 0.01
itsea_b86 0.72 0.25 0.03 0.00
itsea_e5 0.77 0.21 0.02 0.00
itsea_e6 0.81 0.18 0.01 0.00
itsea_e11 0.92 0.07 0.01 0.00
itsea_a7 0.72 0.23 0.05 0.00
itsea_a21 0.58 0.40 0.02 0.01
itsea_a23 0.60 0.35 0.05 0.00
itsea_b31 0.10 0.77 0.13 0.00
itsea_b45 0.77 0.22 0.01 0.00
itsea_b50 0.69 0.27 0.04 0.00
itsea_b53 0.72 0.28 0.00 0.00
itsea_b59 0.42 0.51 0.07 0.00
itsea_b65 0.11 0.73 0.16 0.00
itsea_b66 0.82 0.16 0.02 0.01
itsea_b74 0.77 0.23 0.00 0.02
itsea_b80 0.78 0.17 0.05 0.01
itsea_b85 0.55 0.44 0.01 0.00
d0 %>%
summarize_at(
vars(
age_18mo,
age_lena_d1,
NEG,
duration_total,
duration_analyzed,
ct_first_time,
ct_last_time,
cvc_hour_max,
ctc_hour_max,
awc_hour_max,
ctc_prop,
awc_prop,
sens_r_sfp,
sens_fp,
negmood_fp,
itsea_symptoms,
mom_age,
income_needs,
percent_mother,
cesd_t3,
cesd_t4
),
funs(mean, sd, min, max), na.rm = TRUE
)
`funs()` is deprecated as of dplyr 0.8.0.
Please use a list of either functions or lambdas:
# Simple named list:
list(mean = mean, median = median)
# Auto named with `tibble::lst()`:
tibble::lst(mean, median)
# Using lambdas
list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.
d0 %>%
count(male) %>%
mutate(
per = n / sum(n)
)
d0 %>%
count(itsea_concern_ied) %>%
mutate(
per = n / sum(n)
)
d0 %>%
count(income_needs < 1) %>%
mutate(
per = n / sum(n)
)
d0 %>%
count(mom_latinx) %>%
mutate(
per = n / sum(n)
)
d0 %>%
count(mom_race) %>%
mutate(
per = n / sum(n)
) %>%
arrange(desc(n))
d0 %>%
count(education_txt) %>%
mutate(
per = n / sum(n)
) %>%
arrange(education_txt)
d0 %>%
count(day_type) %>%
mutate(per = n / sum(n))
d0 %>%
count(cesd_t3 >= 16) %>%
mutate(per = n / sum(n))
d0 %>%
count(cesd_t4 >= 16) %>%
mutate(per = n / sum(n))
NA
d0 %>%
count(lena_mornstart, lena_recordfull)
d0 %>%
filter(lena_mornstart == 0) %>%
dplyr::select(
ID,
lena_recordno,
ct_first_time,
ct_last_time,
duration_total,
duration_analyzed,
lena_notes
)
d0 %>%
identify_outliers_hist(x = cvc_hour_max)
d0 %>%
identify_outliers_hist(x = ctc_hour_max)
d0 %>%
identify_outliers_hist(x = ctc_prop)
d0 %>%
identify_outliers_hist(x = awc_hour_max)
d0 %>%
identify_outliers_hist(x = awc_prop)
d0 %>%
identify_outliers_hist(x = sens_r_sfp)
d0 %>%
identify_outliers_hist(x = sens_fp)
d0 %>%
identify_outliers_hist(x = itsea_symptoms)
d0 %>%
identify_outliers_hist(x = cesd_t3)
d0 %>%
identify_outliers_hist(x = cesd_t4)
d0 %>%
ggplot(aes(percent_mother)) +
geom_histogram() +
theme_lena +
labs(
x = "% of the LENA recording day infant with mother"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/percent_mom_hist.png",
width = 7,
height = 5
)
corr_lena <-
d0 %>%
dplyr::select(
`AW consistency` = awc_prop,
`CT consistency` = ctc_prop,
`AW quantity` = awc_hour_max,
`CT quantity` = ctc_hour_max,
`Infant vocalizations` = cvc_hour_max,
`Sensitivity to distress` = sens_r_sfp,
`Sensitivity during play` = sens_fp,
) %>%
correlate(use = "pairwise.complete.obs", method = "pearson") %>%
shave(upper = TRUE) %>%
fashion()
Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'
corr_lena_plot <-
corr_lena %>%
gather(variable, value, -rowname) %>%
mutate(
variable = str_replace_all(variable, "\\.", " "),
rowname = as.character(rowname),
value_chr = as.character(value),
value_chr = case_when(
rowname == "CT consistency" & variable == "AW quantity" ~ ".23",
rowname == "AW quantity" & variable == "CT consistency" ~ NA_character_,
rowname == "Sensitivity to distress" & variable == "Sensitivity during play" ~ ".10",
rowname == "Sensitivity during play" & variable == "Sensitivity to distress" ~ NA_character_,
TRUE ~ value_chr
),
value_num = as.numeric(value_chr)
)
corr_lena_plot %>%
ggplot(aes(x = rowname, y = variable)) +
geom_tile(aes(fill = abs(value_num))) +
geom_text(
aes(label = value_chr),
size = 4
) +
scale_fill_gradient2(
low = "blue",
high = "red",
mid = "white",
na.value = "white",
midpoint = 0,
limit = c(-1, 1),
space = "Lab",
name = "Pearson correlation\ncoefficient"
) +
theme_void() +
theme(
axis.text.x = element_text(
angle = 45,
hjust = 1,
vjust = 1,
size = 12
),
axis.text.y = element_text(
size = 12,
hjust = 1.1
)
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/lena_corr_plot.eps",
dpi = 1000,
height = 5,
width = 7
)
cor.test(d0$cesd_t3, d0$awc_hour_max)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$awc_hour_max
t = 0.70855, df = 98, p-value = 0.4803
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1268048 0.2641061
sample estimates:
cor
0.07139142
cor.test(d0$cesd_t3, d0$awc_prop)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$awc_prop
t = 0.20234, df = 98, p-value = 0.8401
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1766922 0.2159863
sample estimates:
cor
0.02043509
cor.test(d0$cesd_t3, d0$ctc_hour_max)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$ctc_hour_max
t = 0.30729, df = 98, p-value = 0.7593
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1664066 0.2260661
sample estimates:
cor
0.03102564
cor.test(d0$cesd_t3, d0$ctc_prop)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$ctc_prop
t = -0.040919, df = 98, p-value = 0.9674
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2003889 0.1924409
sample estimates:
cor
-0.00413343
cor.test(d0$cesd_t3, d0$cvc_hour_max)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$cvc_hour_max
t = 0.046727, df = 98, p-value = 0.9628
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1918759 0.2009519
sample estimates:
cor
0.00472013
cor.test(d0$cesd_t3, d0$sens_r_sfp)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$sens_r_sfp
t = 0.75743, df = 92, p-value = 0.4507
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1259026 0.2769225
sample estimates:
cor
0.07872259
cor.test(d0$cesd_t3, d0$sens_fp)
Pearson's product-moment correlation
data: d0$cesd_t3 and d0$sens_fp
t = -0.18924, df = 94, p-value = 0.8503
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2191435 0.1816818
sample estimates:
cor
-0.01951496
cor.test(d0$NEG, d0$awc_hour_max)
Pearson's product-moment correlation
data: d0$NEG and d0$awc_hour_max
t = 0.81245, df = 94, p-value = 0.4186
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1189726 0.2793149
sample estimates:
cor
0.08350513
cor.test(d0$NEG, d0$awc_prop)
Pearson's product-moment correlation
data: d0$NEG and d0$awc_prop
t = -0.52757, df = 94, p-value = 0.599
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2520746 0.1477609
sample estimates:
cor
-0.05433463
cor.test(d0$NEG, d0$ctc_hour_max)
Pearson's product-moment correlation
data: d0$NEG and d0$ctc_hour_max
t = 1.4766, df = 94, p-value = 0.1431
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.05147543 0.34076399
sample estimates:
cor
0.1505644
cor.test(d0$NEG, d0$ctc_prop)
Pearson's product-moment correlation
data: d0$NEG and d0$ctc_prop
t = -0.3159, df = 94, p-value = 0.7528
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2315400 0.1690237
sample estimates:
cor
-0.03256579
cor.test(d0$NEG, d0$cvc_hour_max)
Pearson's product-moment correlation
data: d0$NEG and d0$cvc_hour_max
t = 0.24743, df = 94, p-value = 0.8051
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1758738 0.2248476
sample estimates:
cor
0.02551173
cor.test(d0$NEG, d0$sens_r_sfp)
Pearson's product-moment correlation
data: d0$NEG and d0$sens_r_sfp
t = -0.2727, df = 88, p-value = 0.7857
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2347363 0.1791112
sample estimates:
cor
-0.02905776
cor.test(d0$NEG, d0$sens_fp)
Pearson's product-moment correlation
data: d0$NEG and d0$sens_fp
t = -1.1946, df = 90, p-value = 0.2354
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.3215266 0.0819784
sample estimates:
cor
-0.1249368
cor.test(d0$awc_hour_max, d0$sens_fp)
cor.test(d0$awc_hour_max, d0$sens_r_sfp)
cor.test(d0$ctc_hour_max, d0$sens_fp)
cor.test(d0$ctc_hour_max, d0$sens_r_sfp)
cor.test(d0$awc_prop, d0$sens_fp)
cor.test(d0$awc_prop, d0$sens_r_sfp)
cor.test(d0$ctc_prop, d0$sens_fp)
cor.test(d0$ctc_prop, d0$sens_r_sfp)
cor.test(d0$cvc_hour_max, d0$sens_fp)
cor.test(d0$cvc_hour_max, d0$sens_r_sfp)
cor.test(d0$sens_fp, d0$sens_r_sfp)
t.test(d0$awc_hour_max ~ d0$secondlang)
Welch Two Sample t-test
data: d0$awc_hour_max by d0$secondlang
t = 0.28166, df = 94.674, p-value = 0.7788
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-503.7566 670.3237
sample estimates:
mean in group 0 mean in group 1
3549.149 3465.865
t.test(d0$ctc_hour_max ~ d0$secondlang)
Welch Two Sample t-test
data: d0$ctc_hour_max by d0$secondlang
t = 0.47747, df = 91.869, p-value = 0.6342
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-9.311072 15.204689
sample estimates:
mean in group 0 mean in group 1
75.44681 72.50000
t.test(d0$ctc_prop ~ d0$secondlang)
Welch Two Sample t-test
data: d0$ctc_prop by d0$secondlang
t = 0.87779, df = 96.693, p-value = 0.3822
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.02600516 0.06724637
sample estimates:
mean in group 0 mean in group 1
0.5561746 0.5355540
t.test(d0$awc_prop ~ d0$secondlang)
Welch Two Sample t-test
data: d0$awc_prop by d0$secondlang
t = 0.52872, df = 96.086, p-value = 0.5982
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.03152325 0.05441354
sample estimates:
mean in group 0 mean in group 1
0.7422157 0.7307705
t.test(d0$itsea_symptoms ~ d0$secondlang)
Welch Two Sample t-test
data: d0$itsea_symptoms by d0$secondlang
t = -0.34477, df = 96.903, p-value = 0.731
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.2353910 0.1657158
sample estimates:
mean in group 0 mean in group 1
1.053025 1.087863
t.test(d0$awc_hour_max ~ d0$day_type)
Welch Two Sample t-test
data: d0$awc_hour_max by d0$day_type
t = -1.4866, df = 66.225, p-value = 0.1419
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1072.0120 156.9219
sample estimates:
mean in group weekday mean in group weekend
3348.369 3805.914
t.test(d0$ctc_hour_max ~ d0$day_type)
Welch Two Sample t-test
data: d0$ctc_hour_max by d0$day_type
t = 1.0076, df = 69.607, p-value = 0.3171
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-6.267328 19.062933
sample estimates:
mean in group weekday mean in group weekend
76.36923 69.97143
t.test(d0$ctc_prop ~ d0$day_type)
Welch Two Sample t-test
data: d0$ctc_prop by d0$day_type
t = 0.74512, df = 65.551, p-value = 0.4589
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.03151406 0.06903400
sample estimates:
mean in group weekday mean in group weekend
0.5534388 0.5346789
t.test(d0$awc_prop ~ d0$day_type)
Welch Two Sample t-test
data: d0$awc_prop by d0$day_type
t = -1.4395, df = 68.83, p-value = 0.1545
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.07686127 0.01243104
sample estimates:
mean in group weekday mean in group weekend
0.7255763 0.7577915
d0 %>%
dplyr::select(
ctc_hour_max,
awc_hour_max,
ctc_prop,
awc_prop,
sens_fp,
sens_r_sfp
) %>%
gather(measure, value, ctc_hour_max:sens_r_sfp) %>%
mutate(
measure = factor(
measure,
levels = c(
"awc_prop",
"ctc_prop",
"awc_hour_max",
"ctc_hour_max",
"sens_r_sfp",
"sens_fp"
),
labels = c(
"AW consistency",
"CT consistency",
"AW quantity",
"CT quantity",
"Sensitivity to distress",
"Sensivity during play"
)
)
) %>%
ggplot(aes(value)) +
geom_histogram() +
facet_wrap(.~measure, scales = "free", ncol = 2, nrow = 3) +
theme_minimal() +
theme(
strip.text = element_text(size = 18),
panel.grid = element_blank(),
plot.title = element_text(size = 18, hjust = .5),
axis.title = element_text(size = 18),
axis.text = element_text(size = 16),
legend.position = "none"
) +
labs(
x = "Caregiving input value"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/care_hist_ggpanels.png",
width = 7.5,
height = 8
)
d0 %>%
ggplot(aes(itsea_symptoms)) +
geom_histogram() +
scale_x_continuous(breaks = seq.int(0, 3, .25)) +
theme_lena +
labs(
x = "Toddler symptoms of psychopathology"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/itsea_hist.png"
)
Saving 6.33 x 3.91 in image
d0 %>%
dplyr::select(
Activity = itsea_activity,
Aggression = itsea_agress,
Depression = itsea_depp,
Anxiety = itsea_gad,
`Negative emotion` = itsea_negemo
) %>%
gather(itsea_subscale, score, Activity:`Negative emotion`) %>%
ggplot(aes(score)) +
geom_histogram(bins = 12) +
theme_lena +
theme(
strip.text = element_text(size = 16)
) +
facet_wrap(.~itsea_subscale, scales = "free") +
labs(
x = "Subscale score"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/itsea_subscale_hist.png",
width = 10,
height = 7
)
d0 %>%
ggplot(aes(cesd_t3)) +
geom_histogram() +
scale_x_continuous(breaks = seq.int(0, 50, 5)) +
theme_lena +
labs(
x = "Maternal depressive symptoms\n at infant age 6 months"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/cesd6_hist.png"
)
Saving 6.33 x 3.91 in image
d0 %>%
ggplot(aes(cesd_t4)) +
geom_histogram() +
scale_x_continuous(breaks = seq.int(0, 50, 5)) +
theme_lena +
labs(
x = "Maternal depressive symptoms\n at infant age 18 months"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/cesd18_hist.png"
)
Saving 6.33 x 3.91 in image
d0 <-
d0 %>%
mutate(
white = as.factor(
if_else(
mom_race == "White",
1, 0
)
),
male = as.factor(male),
mom_latinx = as.factor(mom_latinx)
)
contrasts(d0$male) <- c(-.5, .5)
contrasts(d0$white) <- c(-.5, .5)
contrasts(d0$mom_latinx) <- c(-.5, .5)
d0_lm <-
d0 %>%
filter(!is.na(NEG))
mod_0 <-
lm(
scale(itsea_symptoms) ~
1,
data = d0_lm
)
mod_cesdt3 <-
lm(
scale(itsea_symptoms) ~
scale(cesd_t3),
data = d0_lm
)
anova(mod_0, mod_cesdt3)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ 1
Model 2: scale(itsea_symptoms) ~ scale(cesd_t3)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 95 95.000
2 94 93.846 1 1.1538 1.1557 0.2851
mod_cesdt4 <-
lm(
scale(itsea_symptoms) ~
scale(cesd_t4),
data = d0_lm
)
anova(mod_0, mod_cesdt4)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ 1
Model 2: scale(itsea_symptoms) ~ scale(cesd_t4)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 95 95.000
2 94 91.346 1 3.6537 3.7598 0.0555 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
mod_neg <-
lm(
scale(itsea_symptoms) ~
scale(cesd_t4) +
scale(NEG),
data = d0_lm
)
anova(mod_cesdt4, mod_neg)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4)
Model 2: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 94 91.346
2 93 81.119 1 10.227 11.725 0.000919 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
model_parameters(mod_neg)
Parameter | Coefficient | SE | 95% CI | t | df | p
-------------------------------------------------------------------------
(Intercept) | 9.39e-17 | 0.10 | [-0.19, 0.19] | 9.85e-16 | 93 | > .999
cesd_t4 | 0.18 | 0.10 | [-0.01, 0.37] | 1.84 | 93 | 0.069
NEG | 0.33 | 0.10 | [ 0.14, 0.52] | 3.42 | 93 | < .001
calc.relimp(mod_neg)
Response variable: scale(itsea_symptoms)
Total response variance: 1
Analysis based on 96 observations
2 Regressors:
scale(cesd_t4) scale(NEG)
Proportion of variance explained by model: 14.61%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(cesd_t4) 0.03476271
scale(NEG) 0.11135116
Average coefficients for different model sizes:
1X 2Xs
scale(cesd_t4) 0.1961112 0.1765674
scale(NEG) 0.3391873 0.3286887
mod_Bage <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(age_lena_d1),
data = d0_lm
)
anova(mod_neg, mod_Bage)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(age_lena_d1)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 79.246 1 1.8733 2.1748 0.1437
mod_Mage <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(mom_age),
data = d0_lm
)
anova(mod_neg, mod_Mage)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(mom_age)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 78.782 1 2.3375 2.7298 0.1019
mod_race <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
white,
data = d0_lm
)
anova(mod_neg, mod_race)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + white
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 80.305 1 0.81468 0.9333 0.3365
mod_ethnicity <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
mom_latinx,
data = d0_lm
)
anova(mod_neg, mod_ethnicity)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + mom_latinx
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 79.287 1 1.8318 2.1255 0.1483
mod_education <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(education),
data = d0_lm
)
anova(mod_neg, mod_education)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(education)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 79.548 1 1.5715 1.8175 0.1809
mod_income_needs <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(income_needs),
data = d0_lm
)
anova(mod_neg, mod_income_needs)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + scale(income_needs)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 80.995 1 0.12377 0.1406 0.7086
mod_sex <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
male,
data = d0_lm
)
anova(mod_neg, mod_sex)
Analysis of Variance Table
Model 1: scale(itsea_symptoms) ~ scale(cesd_t4) + scale(NEG)
Model 2: scale(itsea_symptoms) ~ scale(NEG) + scale(cesd_t4) + male
Res.Df RSS Df Sum of Sq F Pr(>F)
1 93 81.119
2 92 79.115 1 2.004 2.3304 0.1303
mod_awprop <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(awc_prop),
data = d0_lm
)
model_parameters(mod_awprop)
Parameter | Coefficient | SE | 95% CI | t | df | p
--------------------------------------------------------------------------
(Intercept) | 1.22e-16 | 0.09 | [-0.19, 0.19] | 1.31e-15 | 92 | > .999
NEG | 0.32 | 0.09 | [ 0.13, 0.50] | 3.36 | 92 | 0.001
cesd_t4 | 0.20 | 0.09 | [ 0.01, 0.38] | 2.09 | 92 | 0.040
awc_prop | -0.22 | 0.09 | [-0.41, -0.03] | -2.32 | 92 | 0.023
calc.relimp(mod_awprop)
Response variable: scale(itsea_symptoms)
Total response variance: 1
Analysis based on 96 observations
3 Regressors:
scale(NEG) scale(cesd_t4) scale(awc_prop)
Proportion of variance explained by model: 19.33%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(NEG) 0.10717209
scale(cesd_t4) 0.03851538
scale(awc_prop) 0.04762374
Average coefficients for different model sizes:
1X 2Xs 3Xs
scale(NEG) 0.3391873 0.3284913 0.3156267
scale(cesd_t4) 0.1961112 0.1967953 0.1965872
scale(awc_prop) -0.2183248 -0.2189630 -0.2184895
mod_ctprop <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(ctc_prop),
data = d0_lm
)
model_parameters(mod_ctprop)
Parameter | Coefficient | SE | 95% CI | t | df | p
--------------------------------------------------------------------------
(Intercept) | 1.63e-17 | 0.09 | [-0.18, 0.18] | 1.76e-16 | 92 | > .999
NEG | 0.32 | 0.09 | [ 0.13, 0.50] | 3.43 | 92 | < .001
cesd_t4 | 0.20 | 0.09 | [ 0.02, 0.39] | 2.16 | 92 | 0.034
ctc_prop | -0.25 | 0.09 | [-0.43, -0.06] | -2.63 | 92 | 0.010
calc.relimp(mod_ctprop)
Response variable: scale(itsea_symptoms)
Total response variance: 1
Analysis based on 96 observations
3 Regressors:
scale(NEG) scale(cesd_t4) scale(ctc_prop)
Proportion of variance explained by model: 20.6%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(NEG) 0.10841060
scale(cesd_t4) 0.03952569
scale(ctc_prop) 0.05808721
Average coefficients for different model sizes:
1X 2Xs 3Xs
scale(NEG) 0.3391873 0.3302660 0.3191677
scale(cesd_t4) 0.1961112 0.1993232 0.2018538
scale(ctc_prop) -0.2363189 -0.2420650 -0.2461931
mod_awmax<-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(awc_hour_max),
data = d0_lm
)
model_parameters(mod_awmax)
Parameter | Coefficient | SE | 95% CI | t | df | p
--------------------------------------------------------------------------
(Intercept) | 8.29e-17 | 0.10 | [-0.19, 0.19] | 8.67e-16 | 92 | > .999
NEG | 0.33 | 0.10 | [ 0.14, 0.53] | 3.45 | 92 | < .001
cesd_t4 | 0.18 | 0.10 | [-0.01, 0.37] | 1.86 | 92 | 0.066
awc_hour_max | -0.06 | 0.10 | [-0.25, 0.13] | -0.62 | 92 | 0.537
calc.relimp(mod_awmax)
Response variable: scale(itsea_symptoms)
Total response variance: 1
Analysis based on 96 observations
3 Regressors:
scale(NEG) scale(cesd_t4) scale(awc_hour_max)
Proportion of variance explained by model: 14.97%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(NEG) 0.112514996
scale(cesd_t4) 0.035192701
scale(awc_hour_max) 0.001947467
Average coefficients for different model sizes:
1X 2Xs 3Xs
scale(NEG) 0.33918734 0.33605409 0.33350200
scale(cesd_t4) 0.19611119 0.18725146 0.17959080
scale(awc_hour_max) -0.02200461 -0.04182114 -0.05979412
mod_ctmax <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(ctc_hour_max),
data = d0_lm
)
model_parameters(mod_ctmax)
Parameter | Coefficient | SE | 95% CI | t | df | p
---------------------------------------------------------------------------
(Intercept) | 5.46e-17 | 0.09 | [-0.19, 0.19] | 5.86e-16 | 92 | > .999
NEG | 0.36 | 0.09 | [ 0.17, 0.55] | 3.79 | 92 | < .001
cesd_t4 | 0.21 | 0.10 | [ 0.02, 0.40] | 2.23 | 92 | 0.028
ctc_hour_max | -0.22 | 0.10 | [-0.41, -0.03] | -2.28 | 92 | 0.025
calc.relimp(mod_ctmax)
Response variable: scale(itsea_symptoms)
Total response variance: 1
Analysis based on 96 observations
3 Regressors:
scale(NEG) scale(cesd_t4) scale(ctc_hour_max)
Proportion of variance explained by model: 19.19%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(NEG) 0.12030047
scale(cesd_t4) 0.04075023
scale(ctc_hour_max) 0.03087917
Average coefficients for different model sizes:
1X 2Xs 3Xs
scale(NEG) 0.3391873 0.3478042 0.3596103
scale(cesd_t4) 0.1961112 0.2007037 0.2123774
scale(ctc_hour_max) -0.1289444 -0.1758479 -0.2195135
mod_cvmax <-
lm(
scale(itsea_symptoms) ~
scale(NEG) +
scale(cesd_t4) +
scale(cvc_hour_max),
data = d0_lm
)
model_parameters(mod_cvmax)
Parameter | Coefficient | SE | 95% CI | t | df | p
---------------------------------------------------------------------------
(Intercept) | 5.25e-17 | 0.09 | [-0.18, 0.18] | 5.67e-16 | 92 | > .999
NEG | 0.33 | 0.09 | [ 0.15, 0.52] | 3.56 | 92 | < .001
cesd_t4 | 0.22 | 0.09 | [ 0.03, 0.40] | 2.28 | 92 | 0.025
cvc_hour_max | -0.24 | 0.09 | [-0.43, -0.05] | -2.54 | 92 | 0.013
calc.relimp(mod_cvmax)
Response variable: scale(itsea_symptoms)
Total response variance: 1
Analysis based on 96 observations
3 Regressors:
scale(NEG) scale(cesd_t4) scale(cvc_hour_max)
Proportion of variance explained by model: 20.21%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(NEG) 0.11276570
scale(cesd_t4) 0.04185695
scale(cvc_hour_max) 0.04746378
Average coefficients for different model sizes:
1X 2Xs 3Xs
scale(NEG) 0.3391873 0.3365568 0.3325144
scale(cesd_t4) 0.1961112 0.2052806 0.2151002
scale(cvc_hour_max) -0.1965144 -0.2198211 -0.2397697
model_d <-
d0_lm %>%
dplyr::select(
awc_prop,
ctc_prop,
awc_hour_max,
ctc_hour_max,
cvc_hour_max,
itsea_symptoms,
NEG,
cesd_t4
) %>%
as.data.frame()
# each measure is compared to a covariate only model
lmBF_covonly = lmBF(itsea_symptoms ~ NEG + cesd_t4, data = model_d)
lmBF_awcprop = lmBF(itsea_symptoms ~ awc_prop + cesd_t4 + NEG, data = model_d)
lmBF_ctcprop = lmBF(itsea_symptoms ~ ctc_prop + cesd_t4 + NEG, data = model_d)
lmBF_awcmax = lmBF(itsea_symptoms ~ awc_hour_max + cesd_t4 + NEG, data = model_d)
lmBF_ctcmax = lmBF(itsea_symptoms ~ ctc_hour_max + cesd_t4 + NEG, data = model_d)
mBF_cvcmax = lmBF(itsea_symptoms ~ cvc_hour_max + cesd_t4 + NEG, data = model_d)
# BF for AW consistency
lmBF_awcprop / lmBF_covonly
Bayes factor analysis
--------------
[1] awc_prop + cesd_t4 + NEG : 3.008407 ±0%
Against denominator:
itsea_symptoms ~ NEG + cesd_t4
---
Bayes factor type: BFlinearModel, JZS
# BF for CT consistency
lmBF_ctcprop / lmBF_covonly
Bayes factor analysis
--------------
[1] ctc_prop + cesd_t4 + NEG : 5.794294 ±0%
Against denominator:
itsea_symptoms ~ NEG + cesd_t4
---
Bayes factor type: BFlinearModel, JZS
# BF for AW quantity
lmBF_awcmax / lmBF_covonly
Bayes factor analysis
--------------
[1] awc_hour_max + cesd_t4 + NEG : 0.3489947 ±0%
Against denominator:
itsea_symptoms ~ NEG + cesd_t4
---
Bayes factor type: BFlinearModel, JZS
# BF for CT quantity
lmBF_ctcmax / lmBF_covonly
Bayes factor analysis
--------------
[1] ctc_hour_max + cesd_t4 + NEG : 2.803794 ±0%
Against denominator:
itsea_symptoms ~ NEG + cesd_t4
---
Bayes factor type: BFlinearModel, JZS
# BF for infant vocalizations
mBF_cvcmax / lmBF_covonly
Bayes factor analysis
--------------
[1] cvc_hour_max + cesd_t4 + NEG : 4.723113 ±0%
Against denominator:
itsea_symptoms ~ NEG + cesd_t4
---
Bayes factor type: BFlinearModel, JZS
best_cov <- lm(
itsea_symptoms ~
scale(cesd_t4, scale = FALSE) +
scale(NEG, scale = FALSE),
data = d0
)
d0 <-
d0 %>%
add_residuals(best_cov, var = "resid_symptoms")
d0 %>%
ggplot(aes(awc_prop, resid_symptoms)) +
geom_point(size = 4, alpha = 1/2) +
geom_smooth(method = "lm", size = 3, color = "black") +
scale_x_continuous(
breaks = seq.int(0, 1, .10),
labels = scales::percent_format(accuracy = 1)
) +
theme_minimal() +
theme(
strip.text = element_text(size = 17),
panel.grid = element_blank(),
plot.title = element_text(size = 18, hjust = .5),
axis.title = element_text(size = 18),
axis.text = element_text(size = 16),
legend.position = "none"
) +
labs(
x = "Consistency of adult words in daily life",
y = "Toddler symptoms of psychopathology\n(residuals)"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/aw_prop_symptoms.png",
width = 7,
height = 6
)
d0 %>%
ggplot(aes(ctc_prop, resid_symptoms)) +
geom_point(size = 4, alpha = 1/2) +
geom_smooth(method = "lm", size = 3, color = "black") +
scale_x_continuous(
breaks = seq.int(0, 1, .10),
labels = scales::percent_format(accuracy = 1)
) +
theme_minimal() +
theme(
strip.text = element_text(size = 17),
panel.grid = element_blank(),
plot.title = element_text(size = 18, hjust = .5),
axis.title = element_text(size = 18),
axis.text = element_text(size = 16),
legend.position = "none"
) +
labs(
x = "Consistency of conversational turns in daily life",
y = "Toddler symptoms of psychopathology\n(residuals)"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/ctc_prop_symptoms.png",
width = 7,
height = 6
)
d0 %>%
dplyr::select(ctc_prop, awc_prop, resid_symptoms) %>%
rename(
`CT consistency` = ctc_prop,
`AW consistency` = awc_prop
) %>%
gather(measure, value, -resid_symptoms) %>%
ggplot(aes(value, resid_symptoms)) +
geom_point(size = 4, alpha = 1/2) +
geom_smooth(method = "lm", size = 3, color = "black") +
scale_x_continuous(
breaks = seq.int(0, 1, .15),
labels = scales::percent_format(accuracy = 1)
) +
theme_minimal() +
theme(
strip.text = element_text(size = 17),
panel.grid = element_blank(),
plot.title = element_text(size = 18, hjust = .5),
axis.title = element_text(size = 18),
axis.text = element_text(size = 16),
legend.position = "none"
) +
labs(
x = "Consistency of language input in daily life",
y = "Toddler symptoms of psychopathology\n(residuals)"
) +
facet_grid(.~measure, scales = "free")
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/ctc_awc_prop_symptoms.eps",
width = 10,
height = 6,
dpi = 1000
)
d0 %>%
ggplot(aes(cvc_hour_max, resid_symptoms)) +
geom_point(size = 4, alpha = 1/2) +
geom_smooth(method = "lm", size = 3, color = "black") +
scale_x_continuous(
breaks = seq.int(0, 1000, 100)
) +
theme_minimal() +
theme(
strip.text = element_text(size = 17),
panel.grid = element_blank(),
plot.title = element_text(size = 18, hjust = .5),
axis.title = element_text(size = 18),
axis.text = element_text(size = 16),
legend.position = "none"
) +
labs(
x = "Quantity of infant vocalizations in daily life",
y = "Toddler symptoms of psychopathology\n(residuals)"
)
ggsave(
"~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/plots/cv_max_symptoms.png",
width = 7,
height = 6
)
z_score <- function(x) {
diff_mu <- x - mean(x, na.rm = T)
sd <- sd(x, na.rm = T)
diff_mu / sd
}
d0_en <-
d0 %>%
dplyr::select(
itsea_symptoms,
ctc_prop,
awc_prop,
ctc_hour_max,
awc_hour_max,
cvc_hour_max,
cesd_t4,
NEG
) %>%
mutate_at(
vars(itsea_symptoms:NEG),
funs(z_score)
) %>%
na.omit()
predictors <-
d0_en %>%
dplyr::select(-itsea_symptoms) %>%
as.matrix()
symptoms <- z_score(d0_en$itsea_symptoms)
fit_net <- glmnet(predictors, symptoms, family = "gaussian")
plot(fit_net, label = TRUE)
If nfolds is set to the sample size = leave-one-out CV https://sciphy-stats.com/post/2019-01-25-finalizing-glmnet-models/ Running the CV 100 times to minimize simulation error.
# run cross validation
set.seed(123)
lambdaMins <- c()
for (i in 1:100) {
fit_cv <- cv.glmnet(
predictors,
symptoms,
type.measure = "mse",
family = "gaussian",
alpha = .5,
nfolds = 96,
grouped = FALSE
)
lambdaMins <- cbind(lambdaMins, fit_cv$lambda.min)
}
# save lambda corresponding to minimal mse to a tibble
lambdaMins <- as_tibble(lambdaMins)
The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
Using compatibility `.name_repair`.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.
lambdaMins <-
lambdaMins %>%
gather(run, lambdaMin, V1:V100)
# calculate mean lambda value
lambdaFinal <- mean(lambdaMins$lambdaMin)
fit_net_final <- glmnet(predictors, symptoms, family = "gaussian", alpha = .5, lambda = lambdaFinal)
coefs <- coef(fit_net_final)
en_coef_tbl <-
tidy(coef(fit_net_final)) %>%
dplyr::select(-column) %>%
rename(
Variable = row,
`Elastic net estimate` = value
) %>%
mutate(
Variable = recode(
Variable,
"(Intercept)" = "(Intercept)",
"ctc_prop" = "CT consistency",
"awc_prop" = "AW consistency",
"ctc_hour_max" = "CT quantity",
"cvc_hour_max" = "Infant vocalizations",
"cesd_t4" = "Maternal depressive symptoms",
"NEG" = "Infant negative affectivity"
)
) %>%
mutate_at(vars(`Elastic net estimate`) , round, 2)
'tidy.dgCMatrix' is deprecated.
See help("Deprecated")'tidy.dgTMatrix' is deprecated.
See help("Deprecated")
# save model matrix
X <- model.matrix(itsea_symptoms ~ ., d0_en)
# get predicted values
y_hat <- matrix(X, ncol = 8) %*% as.matrix(coefs)
# calculate R-squared
r_sq <- 1 - (sum((y_hat[, 1] - symptoms)^2) / sum((symptoms - mean(symptoms))^2))
for (i in 1000) {
fit_rr_cv <- cv.glmnet(predictors, symptoms, type.measure = "mse", alpha = .5, family = "gaussian")
lambdaMins <- cbind(lambdaMins, fit_rr_cv$lambda.min)
}
sym_aw_ms <-
lm(
scale(itsea_symptoms) ~
scale(awc_prop) +
scale(cesd_t4) +
scale(NEG) +
scale(sens_fp) +
scale(sens_r_sfp),
data = d0
)
calc.relimp(sym_aw_ms)
Response variable: scale(itsea_symptoms)
Total response variance: 1.019533
Analysis based on 89 observations
5 Regressors:
scale(awc_prop) scale(cesd_t4) scale(NEG) scale(sens_fp) scale(sens_r_sfp)
Proportion of variance explained by model: 24.51%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(awc_prop) 0.052732752
scale(cesd_t4) 0.028702483
scale(NEG) 0.134788698
scale(sens_fp) 0.021700415
scale(sens_r_sfp) 0.007145296
Average coefficients for different model sizes:
1X 2Xs 3Xs 4Xs 5Xs
scale(awc_prop) -0.24321526 -0.23740095 -0.23135613 -0.22535320 -0.2196113
scale(cesd_t4) 0.17364596 0.16915053 0.16557711 0.16284424 0.1608910
scale(NEG) 0.40773727 0.39466127 0.38145121 0.36835829 0.3556209
scale(sens_fp) -0.19538884 -0.17194599 -0.14911310 -0.12717659 -0.1064187
scale(sens_r_sfp) 0.08106572 0.08337868 0.08447331 0.08445072 0.0834660
model_parameters(sym_aw_ms)
Parameter | Coefficient | SE | 95% CI | t | df | p
-----------------------------------------------------------------------
(Intercept) | -0.04 | 0.10 | [-0.23, 0.15] | -0.44 | 83 | 0.662
awc_prop | -0.22 | 0.10 | [-0.41, -0.03] | -2.26 | 83 | 0.026
cesd_t4 | 0.16 | 0.09 | [-0.03, 0.35] | 1.69 | 83 | 0.094
NEG | 0.36 | 0.10 | [ 0.16, 0.56] | 3.54 | 83 | < .001
sens_fp | -0.11 | 0.10 | [-0.31, 0.10] | -1.04 | 83 | 0.300
sens_r_sfp | 0.08 | 0.09 | [-0.11, 0.27] | 0.88 | 83 | 0.382
sym_ct_ms <-
lm(
scale(itsea_symptoms) ~
scale(ctc_prop) +
scale(cesd_t4) +
scale(NEG) +
scale(sens_fp) +
scale(sens_r_sfp),
data = d0
)
calc.relimp(sym_ct_ms)
Response variable: scale(itsea_symptoms)
Total response variance: 1.019533
Analysis based on 89 observations
5 Regressors:
scale(ctc_prop) scale(cesd_t4) scale(NEG) scale(sens_fp) scale(sens_r_sfp)
Proportion of variance explained by model: 25.17%
Metrics are not normalized (rela=FALSE).
Relative importance metrics:
lmg
scale(ctc_prop) 0.055000215
scale(cesd_t4) 0.029029047
scale(NEG) 0.140378673
scale(sens_fp) 0.020233226
scale(sens_r_sfp) 0.007072563
Average coefficients for different model sizes:
1X 2Xs 3Xs 4Xs 5Xs
scale(ctc_prop) -0.23739588 -0.23580976 -0.23431447 -0.2331960 -0.23268930
scale(cesd_t4) 0.17364596 0.16932768 0.16620912 0.1643043 0.16364810
scale(NEG) 0.40773727 0.39784264 0.38836046 0.3795487 0.37165793
scale(sens_fp) -0.19538884 -0.16943354 -0.14329989 -0.1171948 -0.09132713
scale(sens_r_sfp) 0.08106572 0.08349367 0.08442111 0.0838533 0.08183789
model_parameters(sym_ct_ms)
Parameter | Coefficient | SE | 95% CI | t | df | p
-----------------------------------------------------------------------
(Intercept) | -0.05 | 0.10 | [-0.24, 0.14] | -0.56 | 83 | 0.579
ctc_prop | -0.23 | 0.10 | [-0.42, -0.04] | -2.43 | 83 | 0.017
cesd_t4 | 0.16 | 0.09 | [-0.02, 0.35] | 1.73 | 83 | 0.087
NEG | 0.37 | 0.10 | [ 0.17, 0.57] | 3.73 | 83 | < .001
sens_fp | -0.09 | 0.10 | [-0.29, 0.11] | -0.89 | 83 | 0.374
sens_r_sfp | 0.08 | 0.09 | [-0.11, 0.27] | 0.87 | 83 | 0.389
cor.test(d0$sens_fp, d0$itsea_symptoms)
Pearson's product-moment correlation
data: d0$sens_fp and d0$itsea_symptoms
t = -1.8456, df = 94, p-value = 0.0681
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.37348462 0.01401018
sample estimates:
cor
-0.187001
cor.test(d0$sens_r_sfp, d0$itsea_symptoms)
Pearson's product-moment correlation
data: d0$sens_r_sfp and d0$itsea_symptoms
t = 0.79113, df = 92, p-value = 0.4309
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1224545 0.2801529
sample estimates:
cor
0.08220199
d0_follow <-
d0 %>%
filter(
lena_caregivers___1 != 0 |
lena_caregivers___2 != 0 |
lena_caregivers___3 != 0 |
lena_caregivers___4 != 0 |
lena_caregivers___5 != 0 |
lena_caregivers___6 != 0 |
lena_caregivers___7 != 0 |
lena_caregivers___8 != 0 |
lena_caregivers___9 != 0 |
lena_caregivers___10 != 0
) %>%
rename(
mom_present = lena_caregivers___1,
dad_present = lena_caregivers___2,
daycare_present = lena_caregivers___3,
nanny_present = lena_caregivers___4,
sibling_present = lena_caregivers___5,
cousin_present = lena_caregivers___6,
grandparent_present = lena_caregivers___7,
friends_present = lena_caregivers___8,
relative_present = lena_caregivers___9, #any other relative
other_present = lena_caregivers___10
)
d0_follow %>%
count(!is.na(percent_mother))
d0_follow %>%
summarise_at(
vars(percent_mother),
funs(mean, median, sd, min, max), na.rm = TRUE
)
cor.test(d0_follow$percent_mother, d0_follow$ctc_prop, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$percent_mother and d0_follow$ctc_prop
S = 81475, p-value = 0.0774
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.2065836
cor.test(d0_follow$percent_mother, d0_follow$awc_prop, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$percent_mother and d0_follow$awc_prop
S = 72960, p-value = 0.4954
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.08049238
cor.test(d0_follow$percent_mother, d0_follow$ctc_hour_max, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$percent_mother and d0_follow$ctc_hour_max
S = 85137, p-value = 0.0248
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.2608283
cor.test(d0_follow$percent_mother, d0_follow$awc_hour_max, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$percent_mother and d0_follow$awc_hour_max
S = 76045, p-value = 0.2841
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.1261794
cor.test(d0_follow$percent_mother, d0_follow$cvc_hour_max, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$percent_mother and d0_follow$cvc_hour_max
S = 78635, p-value = 0.1613
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.164533
cor.test(d0_follow$percent_mother, d0_follow$itsea_symptoms, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$percent_mother and d0_follow$itsea_symptoms
S = 61876, p-value = 0.4786
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.0836531
d0_follow %>%
count(mom_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(dad_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(mom_present, dad_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(daycare_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(nanny_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(grandparent_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(friends_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(relative_present) %>%
mutate(
per = n / sum(n)
)
d0_follow %>%
count(
mom_present,
dad_present,
daycare_present,
nanny_present,
grandparent_present
) %>%
mutate(
per = n / sum(n)
) %>%
arrange(n)
NA
d0_follow %>%
count(percent_mother) %>%
arrange(desc(n))
d0_follow <-
d0_follow %>%
mutate(
n_caregivers =
mom_present +
dad_present +
daycare_present +
grandparent_present +
nanny_present +
friends_present +
relative_present
)
cor.test(d0_follow$ctc_prop, d0_follow$n_caregivers, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$ctc_prop and d0_follow$n_caregivers
S = 58928, p-value = 0.09238
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.194426
cor.test(d0_follow$awc_prop, d0_follow$n_caregivers, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$awc_prop and d0_follow$n_caregivers
S = 61071, p-value = 0.154
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.1651252
cor.test(d0_follow$ctc_hour_max, d0_follow$n_caregivers, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$ctc_hour_max and d0_follow$n_caregivers
S = 69248, p-value = 0.6472
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.05334703
cor.test(d0_follow$awc_hour_max, d0_follow$n_caregivers, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$awc_hour_max and d0_follow$n_caregivers
S = 46296, p-value = 0.001106
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.3671023
cor.test(d0_follow$cvc_hour_max, d0_follow$n_caregivers, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$cvc_hour_max and d0_follow$n_caregivers
S = 72778, p-value = 0.9653
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.005081446
cor.test(d0_follow$itsea_symptoms, d0_follow$n_caregivers, method = "spearman")
Cannot compute exact p-value with ties
Spearman's rank correlation rho
data: d0_follow$itsea_symptoms and d0_follow$n_caregivers
S = 80673, p-value = 0.3767
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.1028426
d0 %>%
dplyr::select(
ctc_prop,
awc_prop,
ctc_hour_max,
awc_hour_max,
duration_analyzed,
duration_total,
ct_first_time,
ct_last_time
) %>%
correlate(method = "spearman") %>%
fashion()
Correlation method: 'spearman'
Missing treated using: 'pairwise.complete.obs'
cor.test(d0$ctc_prop, d0$duration_analyzed, method = "pearson")
Pearson's product-moment correlation
data: d0$ctc_prop and d0$duration_analyzed
t = -0.8847, df = 98, p-value = 0.3785
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2805269 0.1093159
sample estimates:
cor
-0.08901346
cor.test(d0$awc_prop, d0$duration_analyzed, method = "pearson")
Pearson's product-moment correlation
data: d0$awc_prop and d0$duration_analyzed
t = -2.4278, df = 98, p-value = 0.01702
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.41518015 -0.04381767
sample estimates:
cor
-0.2381858
cor.test(d0$ctc_hour_max, d0$duration_analyzed, method = "pearson")
Pearson's product-moment correlation
data: d0$ctc_hour_max and d0$duration_analyzed
t = 1.509, df = 98, p-value = 0.1345
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.04712287 0.33712968
sample estimates:
cor
0.15069
cor.test(d0$awc_hour_max, d0$duration_analyzed, method = "pearson")
Pearson's product-moment correlation
data: d0$awc_hour_max and d0$duration_analyzed
t = 1.3002, df = 98, p-value = 0.1966
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.06793023 0.31849685
sample estimates:
cor
0.1302255
cor.test(d0$ctc_prop, d0$ct_first_time, method = "pearson")
Pearson's product-moment correlation
data: d0$ctc_prop and d0$ct_first_time
t = -1.7196, df = 98, p-value = 0.08865
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.35561061 0.02615038
sample estimates:
cor
-0.1711468
cor.test(d0$awc_prop, d0$ct_first_time, method = "pearson")
Pearson's product-moment correlation
data: d0$awc_prop and d0$ct_first_time
t = -0.20145, df = 98, p-value = 0.8408
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2159007 0.1767792
sample estimates:
cor
-0.02034532
cor.test(d0$ctc_hour_max, d0$ct_first_time, method = "pearson")
Pearson's product-moment correlation
data: d0$ctc_hour_max and d0$ct_first_time
t = -2.1017, df = 98, p-value = 0.03814
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.38825834 -0.01173866
sample estimates:
cor
-0.2076779
cor.test(d0$awc_hour_max, d0$ct_first_time, method = "pearson")
Pearson's product-moment correlation
data: d0$awc_hour_max and d0$ct_first_time
t = -2.0029, df = 98, p-value = 0.04795
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.37992144 -0.00195945
sample estimates:
cor
-0.1983012
cor.test(d0$ctc_prop, d0$ct_last_time, method = "pearson")
Pearson's product-moment correlation
data: d0$ctc_prop and d0$ct_last_time
t = -2.6758, df = 98, p-value = 0.008739
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.43505086 -0.06799646
sample estimates:
cor
-0.2609297
cor.test(d0$awc_prop, d0$ct_last_time, method = "pearson")
Pearson's product-moment correlation
data: d0$awc_prop and d0$ct_last_time
t = -2.4694, df = 98, p-value = 0.01526
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.41854994 -0.04788765
sample estimates:
cor
-0.2420292
cor.test(d0$ctc_hour_max, d0$ct_last_time, method = "pearson")
Pearson's product-moment correlation
data: d0$ctc_hour_max and d0$ct_last_time
t = -0.92389, df = 98, p-value = 0.3578
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.2841549 0.1054191
sample estimates:
cor
-0.09292313
cor.test(d0$awc_hour_max, d0$ct_last_time, method = "pearson")
Pearson's product-moment correlation
data: d0$awc_hour_max and d0$ct_last_time
t = -1.3196, df = 98, p-value = 0.19
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.32023976 0.06599806
sample estimates:
cor
-0.1321329
cor.test(d0$ct_last_time, d0$ct_first_time, method = "pearson")
Pearson's product-moment correlation
data: d0$ct_last_time and d0$ct_first_time
t = 6.7365, df = 98, p-value = 0.000000001118
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4116546 0.6834769
sample estimates:
cor
0.5625842
d0 %>%
ggplot(aes(duration_analyzed, awc_prop)) +
geom_point() +
geom_smooth(method = "lm")
d0 %>%
ggplot(aes(ct_first_time, awc_hour_max)) +
geom_point() +
geom_smooth(method = "lm")
d0 %>%
ggplot(aes(ct_first_time, ctc_hour_max)) +
geom_point() +
geom_smooth(method = "lm")
d0 %>%
ggplot(aes(ct_last_time, awc_prop)) +
geom_point() +
geom_smooth(method = "lm")
d0 %>%
ggplot(aes(ct_last_time, ctc_prop)) +
geom_point() +
geom_smooth(method = "lm")
cor.test(d0$ctc_rate, d0$itsea_symptoms)
Pearson's product-moment correlation
data: d0$ctc_rate and d0$itsea_symptoms
t = -1.8434, df = 98, p-value = 0.06829
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.36631497 0.01384722
sample estimates:
cor
-0.1830688
cor.test(d0$awc_rate, d0$itsea_symptoms)
Pearson's product-moment correlation
data: d0$awc_rate and d0$itsea_symptoms
t = -1.5658, df = 98, p-value = 0.1206
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.34214889 0.04146029
sample estimates:
cor
-0.1562301
#write_csv(d0, "~/Box/lucy_king_files/BABIES/lena_symptoms/manuscript/dev_science/king_analyzed_data_20200528.csv")